/***************************************************************************
 *
 *   Copyright (C) 2005-2008 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#include "Pulsar/PolnProfileStats.h"
#include "Pulsar/PolnProfile.h"

#include "Pulsar/OnPulseThreshold.h"
#include "Pulsar/GaussianBaseline.h"

using namespace std;

// #define _DEBUG 1

//! Default constructor
Pulsar::PolnProfileStats::PolnProfileStats (const PolnProfile* _profile) try
{
  avoid_zero_determinant = false;
  regions_set = false;
  stats = new ProfileStats;

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::ctor this=" << this
       << " calling select_profile" << endl;
#endif

  select_profile (_profile);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::ctor";
   }

//! Destructor
Pulsar::PolnProfileStats::~PolnProfileStats()
{
}

void Pulsar::PolnProfileStats::set_avoid_zero_determinant (bool flag)
{
  avoid_zero_determinant = flag;
}

//! Set the PolnProfile from which statistics will be derived
void Pulsar::PolnProfileStats::set_profile (const PolnProfile* _profile) try
{
  profile = _profile;

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::set_profile this=" << this
       << " calling build" << endl;
#endif
  build ();
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::set_profile";
   }

//! Set the PolnProfile from which baseline and onpulse mask will be selected
/*! It is assumed that all subsequent PolnProfile instances passed to
  set_profile will have the same phase as set_profile */
void Pulsar::PolnProfileStats::select_profile (const PolnProfile* _profile) try
{
#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::select_profile this=" << this
       << " PolnProfile=" << _profile << endl;
#endif

  profile = _profile;
  regions_set = false;
  build ();
  if (_profile)
    regions_set = true;
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::select_profile (PolnProfile*)";
   }

void Pulsar::PolnProfileStats::select_profile (const Profile* total) try
{
#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::select_profile this=" << this
       << " Profile=" << total << endl;
#endif

  profile = 0;
  stats->select_profile( total );
  regions_set = true;
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::select_profile (Profile*)";
   }

//! Set the on-pulse and baseline regions
void Pulsar::PolnProfileStats::set_regions (const PhaseWeight& on,
					    const PhaseWeight& off) try
{
  stats->set_regions (on, off);
  regions_set = true;
  build ();
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::set_regions";
   }

//! Set the on-pulse and baseline regions
void Pulsar::PolnProfileStats::get_regions (PhaseWeight& on, 
					    PhaseWeight& off) const try
{
  stats->get_regions (on, off);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_regions";
   }

//! Returns the total flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_intensity () const try
{
  stats->set_profile( profile->get_Profile(0) );
  return stats->get_total ();
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_intensity";
   }

//! Returns the total polarized flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_polarized () const try
{
  Profile polarized;
  profile->get_polarized (&polarized);

  stats->set_profile (&polarized);

  return stats->get_total (false);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_polarized";
   }

//! Returns the total polarized flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_polarized_squared () const try
{
  Profile polarized;
  profile->get_polarized_squared (&polarized);

  stats->set_profile (&polarized);

  return stats->get_total (false);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_polarized_squared";
   }

//! Returns the total linearly polarized flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_linear () const try
{
  Profile linear;
  profile->get_linear (&linear);

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::get_total_linear this=" << this << 
    " calling ProfileStats::set_profile regions_set=" << regions_set << endl;
#endif

  stats->set_profile (&linear);

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::get_total_linear this=" << this <<
    " calling ProfileStats::get_total" << endl;
#endif

  return stats->get_total (false);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_linear";
   }

//! Returns the total linearly polarized flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_linear_squared () const try
{
  Profile linear;
  profile->get_linear_squared (&linear);

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::get_total_linear_squared  this=" << this << 
    " calling ProfileStats::set_profile regions_set=" << regions_set << endl;
#endif

  stats->set_profile (&linear);

#if _DEBUG
  cerr << "Pulsar::PolnProfileStats::get_total_linear_squared  this=" << this <<
    " calling ProfileStats::get_total" << endl;
#endif

  return stats->get_total (false);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_linear_squared ";
   }

//! Returns the total circularly polarized flux of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_circular () const try
{
  stats->set_profile (profile->get_Profile(3));
  return stats->get_total ();
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_circular";
   }

//! Returns the total absolute value of circularly polarized flux
Estimate<double> Pulsar::PolnProfileStats::get_total_abs_circular () const try
{
  Profile circular;
  profile->get_circular (&circular);

  stats->set_profile (&circular);

  return stats->get_total (false);
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::total_abs_circular";
   }

//! Get the Stokes parameters of the specified phase bin
Stokes< Estimate<double> >
Pulsar::PolnProfileStats::get_stokes (unsigned ibin) const try
{
  Stokes< Estimate<double> > result = profile->get_Stokes (ibin);

  for (unsigned ipol=0; ipol < 4; ipol++)
    result[ipol].set_variance( get_baseline_variance(ipol).get_value() );

  return result;
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_stokes";
   }

//! Returns the total determinant of the on-pulse phase bins
Estimate<double> Pulsar::PolnProfileStats::get_total_determinant () const try
{
  Estimate<double> total_det;

  for (unsigned ibin=0; ibin < profile->get_nbin(); ibin++)
    if (stats->get_onpulse(ibin))
      total_det += invariant( get_stokes(ibin) );

#if 0
  cerr << "1: " << total_det << endl;

  Profile invint;
  profile->invint (&invint, false);

  stats->set_profile (&invint);

  total_det = stats->get_total (false);

  cerr << "2: " << total_det << endl;
#endif

  return total_det;
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_total_determinant";
   }

Estimate<double>
Pulsar::PolnProfileStats::get_baseline_variance (unsigned ipol) const try
{
  if (baseline_variance[ipol].get_value() == 0)
  {
    stats->set_profile( profile->get_Profile(ipol) );
    baseline_variance[ipol] = stats->get_baseline_variance();
#ifdef _DEBUG
    cerr << "Pulsar::PolnProfileStats::get_baseline_variance ipol="
         << ipol << " var=" << baseline_variance[ipol] << endl;
#endif
  }
  return baseline_variance[ipol];
}
 catch (Error& error)
   {
     throw error += "Pulsar::PolnProfileStats::get_baseline_variance";
   }

void Pulsar::PolnProfileStats::build () try
{
  baseline_variance = Stokes< Estimate<double> > ();

  if (!profile)
    return;

  if (profile->get_state() != Signal::Stokes)
    throw Error (InvalidParam, "Pulsar::PolnProfileStats::build",
		 "input PolnProfile is not in the Stokes state");

  if (!regions_set)
  {
#ifdef _DEBUG
    cerr << "Pulsar::PolnProfileStats::build this=" << this << 
      " call ProfileStats::select_profile"
      " ProfileStats*=" << stats.ptr() << endl;
#endif

    stats->select_profile( profile->get_Profile(0) );

    if (avoid_zero_determinant)
    {
      Profile invint;
      profile->invint (&invint);

      // roughly speaking ...
      double invint_variance = get_baseline_variance(0).val * 4;
      double threshold = 3.0 * sqrt (invint_variance);

#ifdef _DEBUG
      cerr << "before avoid " << stats->get_onpulse_nbin() << endl;
#endif
      stats->deselect_onpulse (&invint, threshold);
#ifdef _DEBUG
      cerr << "after avoid " << stats->get_onpulse_nbin() << endl;
#endif
    }
  }
}
catch (Error& error)
{
  throw error += "Pulsar::PolnProfileStats::build";
}

//! Returns the vairance of the linearly polarized flux
double Pulsar::PolnProfileStats::get_linear_variance () const
{
  return profile->get_linear_variance();
}

#include "Pulsar/PolnProfileStatsInterface.h"

//! Return a text interface that can be used to configure this instance
TextInterface::Parser* Pulsar::PolnProfileStats::get_interface ()
{
  return new Interface (this);
}


